1 Configuración

source(file.path("code","include.R"), encoding = "UTF-8")
dirs$data <- "datosTFG"
dflt$ini <- as.Date("2019-01-01")
dflt$fin <- as.Date("2021-05-17")

2 Introducción

Usamos 2 servidores. Uno principal (juno) y otro para los datos de riego (dibulibu). Veamos lo que hay en juno:

s <- Verdor$conectarServVerdor()
status <- ssh::ssh_exec_wait(s, "ls RiegoInteligente")
## code
## colorRamps
## datos
## disco.sh
## paramsTS
## saves

Tenemos ficheros de datos (datos) y de código (code). Además…

  • hw.sh: devuelve un listado de características del servidor
  • disco.sh: para ver el espacio de disco ocupado
  • colorRamps, paramsTS contienen parámetros de algunos plots.
status <- ssh::ssh_exec_wait(s, "ls RiegoInteligente/datos")
## amb
## geo
## int
## jardines_riego.geojson
## jardines_riego_NDVI_scl_7_8_9_Mean.geojson
## riego
## verdor
status <- ssh::ssh_exec_wait(s, "ls RiegoInteligente/datos/geo")
## Acometidas.csv
## CONTRATOS_AGUA_PARQUES_Y_JARDINES.csv
## Contratos de alta en AQUACIS parques y jardines.xlsx
## IMIDA_estaciones.csv
## jardines.geojson
## Zonas Verdes Agua Potable
## zonasVerdesAP.geojson
status <- ssh::ssh_exec_wait(s, "ls RiegoInteligente/datos/amb")
## et0
## met
status <- ssh::ssh_exec_wait(s, "ls RiegoInteligente/datos/riego")
## agua
## consumos
## consumos.csv
## contadores.csv
## lecturas
## lecturas.csv
## processed.txt
status <- ssh::ssh_exec_wait(s, "ls RiegoInteligente/datos/verdor")
## images
## images.txt
## summaries
## summaries_NDVI_NA_Mean.csv
## summaries_NDVI_scl_7_8_9_Mean.csv

Los datos tienen la estructura descrita en README.Rmd.

status <- ssh::ssh_exec_wait(s, "ls RiegoInteligente/code")
## DatAmb.R
## DatGeo.R
## DatInt.R
## include.R
## logERR.txt
## logOUT.txt
## Riego.R
## run_juno.sh
## run.R
## run.sh
## sen2r
## times.txt
## Verdor.R
ssh::ssh_disconnect(s)

Lo mismo podemos decir de code.

Por otro lado, en dibulibu, está todo lo relacionado con los datos de riego pero también los datos geográficos de jardines, ya que se necesitan los id de los contratos y las áreas de las zonas verdes:

sr <- Riego$conectarServRiego()
status <- ssh::ssh_exec_wait(sr, "ls RiegoInteligente")
## code
## datos
status <- ssh::ssh_exec_wait(sr, "ls RiegoInteligente/datos")
## geo
## riego
status <- ssh::ssh_exec_wait(sr, "ls RiegoInteligente/datos/geo")
## jardines.geojson
status <- ssh::ssh_exec_wait(sr, "ls RiegoInteligente/datos/riego")
## consumos
## consumos.csv
## contadores.csv
## lecturas
## lecturas.csv
## processed.txt
status <- ssh::ssh_exec_wait(sr, "ls RiegoInteligente/code")
## DatGeo.R
## include.R
## logERR.txt
## logOUT.txt
## Riego.R
## run_dibulibu.sh
## run.R
ssh::ssh_disconnect(sr)

3 Datos geográficos (DatGeo.R)

3.1 Datos geográficos de estaciones meteorológicas

Aquí, es de donde se descargó un CSV (paths$geo("estaciones")) con las ubicaciones de las estaciones meteorológicas.

¡IMPORTANTE! La variable dflt$eemm es la que guarda los identificadores de las estaciones que se tienen en cuenta. Si se consiguen datos de más estaciones hay que modificarla para que se tengan en cuenta. Además, hay que asegurarse de que se disponen de los datos de ubicación de dichas estaciones (función DatGeo$leerEstMet).

dflt$eemm
## [1] "CA42" "CA91" "MO12" "MU21" "MU62"
paths$geo("estaciones")
## [1] "./datosTFG/geo/IMIDA_estaciones.csv"
de <- DatGeo$leerEstMet(dflt$eemm)
str(de)
## Classes 'data.table' and 'data.frame':   5 obs. of  8 variables:
##  $ CODEST       : chr  "CA42" "CA91" "MO12" "MU21" ...
##  $ RED          : chr  "SIAR" "SIAM" "SIAR" "SIAR" ...
##  $ PARAJE       : chr  "Balsapintada" "Campillo Abajo" "Pilica" "Los Alamos" ...
##  $ MUNICIPIO    : chr  "Fuente Alamo" "Fuente Alamo" "Torres de Cotillas" "Beniel" ...
##  $ FECHA DE ALTA: chr  "01/09/99" "01/01/96" "27/09/99" "27/09/00" ...
##  $ LATITUD      : num  37.7 37.7 38 38 37.9
##  $ LONGITUD     : num  -1.13 -1.24 -1.3 -1 -1.13
##  $ ALTITUD      : int  138 175 161 27 56
##  - attr(*, ".internal.selfref")=<externalptr>

3.2 Datos geográficos de jardines

3.2.1 Fichero 2: Ubicaciones de los parques

3.2.1.1 Fichero 2.2: Contornos de las zonas verdes zonasVerdesAP.shp (el importante)

paths$geo("zonas")
## [1] "./datosTFG/geo/Zonas Verdes Agua Potable"
zonasTodas <- DatGeo$leerZonasVerdes()
  • Este fichero contiene varios datos de 1713 zonas verdes de agua potable asociadas a 362 acometidas, como el número de polígonos (Poligonos), el área (Area), el perímetro (Perimetro), etc.
  • El identificador de las acometidas se ha renombrado de CAPA_CAD a Acometida. El identificador señalado NO es una clave primaria (hay acometidas con varias zonas verdes).
str(zonasTodas)
## Classes 'sf' and 'data.frame':   1713 obs. of  7 variables:
##  $ Acometida  : chr  "2213606" "2213606" "2213606" "2213937" ...
##  $ Observacion: chr  NA NA NA NA ...
##  $ Area       : num  1624.3 1354.6 161.7 188.3 58.6 ...
##  $ Perimetro  : num  173.6 163.4 55.6 91.2 37.6 ...
##  $ Geometry   :sfc_MULTIPOLYGON of length 1713; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:62, 1:2] 662155 662155 662156 662153 662147 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  $ Poligonos  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Nodos      : int  61 68 16 16 5 22 33 10 6 11 ...
##  - attr(*, "sf_column")= chr "Geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:6] "Acometida" "Observacion" "Area" "Perimetro" ...

Resumimos el fichero anterior, guardando, para cada acometida:

  • su identificador (Acometida),
  • el número de zonas verdes (ZonasVerdes),
  • el número total de polígonos (Poligonos),
  • el área total (Area),
  • el área del polígono más grande (MasGrande),
  • el área del polígono más pequeño (MasPeque),
  • el perímetro total (Perimetro),
  • el perímetro del polígono más largo (MasLargo),
  • el perímetro del polígono más corto (MasCorto).
paths$geo("zonasR")
## [1] "./datosTFG/geo/zonasVerdesAP.geojson"
zonas <- DatGeo$leerResumenZonasVerdes()
str(zonas)
## Classes 'sf' and 'data.frame':   362 obs. of  13 variables:
##  $ Acometida  : chr  "2111920" "2120205" "2120628" "2122863" ...
##  $ ZonasVerdes: int  4 1 1 4 6 1 3 1 2 2 ...
##  $ Poligonos  : int  4 1 1 4 6 1 3 1 2 2 ...
##  $ MasGrande  : num  115.1 77.4 276.5 209.9 81.7 ...
##  $ MasPeque   : num  21.6 77.4 276.5 195.5 21.7 ...
##  $ Area       : num  273.3 77.4 276.5 810.4 271.5 ...
##  $ MasLargo   : num  45.5 43.6 101 60.5 40 ...
##  $ MasCorto   : num  21.2 43.6 101 56.9 18.8 ...
##  $ Perimetro  : num  133.5 43.6 101 235.3 177.7 ...
##  $ MaxNodos   : int  7 5 13 9 9 10 9 5 9 9 ...
##  $ MinNodos   : int  4 5 13 9 5 10 4 5 6 5 ...
##  $ Nodos      : int  22 5 13 36 36 10 18 5 15 14 ...
##  $ Geometry   :sfc_MULTIPOLYGON of length 362; first list element: List of 4
##   ..$ :List of 1
##   .. ..$ : num [1:4, 1:2] 664175 664170 664177 664175 4204834 ...
##   ..$ :List of 1
##   .. ..$ : num [1:7, 1:2] 664177 664173 664176 664181 664183 ...
##   ..$ :List of 1
##   .. ..$ : num [1:4, 1:2] 664188 664193 664185 664188 4204783 ...
##   ..$ :List of 1
##   .. ..$ : num [1:7, 1:2] 664182 664179 664182 664186 664189 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "Geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:12] "Acometida" "ZonasVerdes" "Poligonos" "MasGrande" ...

3.2.1.2 Fichero 2.1: Coordenadas de las acometidas Acometidas.csv

Una acometida se define como la instalación comprendida entre la tubería de la red de distribución (donde hay una válvula) y la instalación interior del edificio (la batería de contadores).

paths$geo("acometidas")
## [1] "./datosTFG/geo/Acometidas.csv"
acometidas <- DatGeo$leerAcometidas(acometidas = zonas$Acometida)
  • Este fichero contiene varios datos asociados a 336 acometidas de contratos, como el diámetro, la longitud, la profundidad, etc.
  • El identificador de las acometidas se ha renombrado de AC_ID a Acometida.
  • Nos dan datos de ubicación precisos: las coordenadas UTM de las acometidas (Abcisa, Norte).

Vemos que hay 336 acometidas, frente a las 362 que había en zonas, por lo que desconocemos la ubicación exacta de 26 acometidas.

str(acometidas)
## Classes 'data.table' and 'data.frame':   336 obs. of  3 variables:
##  $ Acometida: chr  "2213606" "2191239" "2184773" "2183939" ...
##  $ Abcisa   : num  662181 669069 661004 663640 662653 ...
##  $ Norte    : num  4202444 4177798 4199396 4209313 4208150 ...
##  - attr(*, ".internal.selfref")=<externalptr>

3.2.2 Fichero 1: Contratos de agua de parques y jardines

Tenemos varias versiones de este fichero.

3.2.2.1 Fichero 1.2: última versión Contratos de alta en AQUACIS parques y jardines.xlsx (el importante)

paths$geo("relaciones")
## [1] "./datosTFG/geo/Contratos de alta en AQUACIS parques y jardines.xlsx"
relaciones <- DatGeo$leerRelaciones(acometidas = zonas$Acometida)
  • Este fichero contiene 330 contratos de abastecimiento de agua con EMUASA por parte del ayuntamiento de Murcia, principalmente.
  • Cada contrato tiene un identificador único (Contrato).
  • El contrato depende del tipo de instalación (parque, jardín, …) sobre la que se hace (Agrupacion) y de para qué se usa en particular el punto de servicio (TipoPuntoServicios).
  • El consumo es medido por contadores, que tienen un identificador (Contador). También contiene otros datos del contador: el calibre (Calibre), el modelo (ModeloContador), el tipo (TipoContador), cuándo se instaló (FechaInstContador), además de si funciona con un sistema de telelectura (Telectura).
  • En este caso, además del identificador del contrato, aparece el de la acometida (Acometida).
  • También nos dan más datos de ubicación (Direccion, ComplementoDireccion, Localidad, Barrio).

Lo interesante de esta tabla es que relaciona un parque (Contrato) con la acometida (Acometida) y el contador (Contador). La acometida NO es una clave primaria, pues hay 1 acometida que tiene asociada 2 contratos (dejamos de tener en cuenta este parque, ya no aparece), tampoco lo es el contador, pues un contrato no tiene, pero el contrato sí que es clave primaria.

Vemos que hay 330 acometidas, frente a las 362 que había en zonas, por lo que desconocemos el contrato de 362-330=32 acometidas.

str(relaciones)
## Classes 'data.table' and 'data.frame':   330 obs. of  15 variables:
##  $ Contrato            : chr  "6351306" "6351960" "6353617" "6353635" ...
##  $ Acometida           : chr  "2185186" "2183895" "2183457" "2183384" ...
##  $ Actividad           : Factor w/ 4 levels "ACTIVIDAD GENERAL",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Ext/Int             : Factor w/ 2 levels "E","I": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Prop/Inqui          : Factor w/ 3 levels "C","I","P": 1 1 1 3 3 1 1 1 1 1 ...
##  $ Direccion           : chr  "FEDERICO BALART SN C-1-JA" "VICENTE ALEIXANDRE SN C-1-JA" "DEL PILAR SN C-1-JA" "PL.IGLESIA SN C-1-JA" ...
##  $ ComplementoDireccion: chr  NA NA NA NA ...
##  $ Localidad           : Factor w/ 49 levels "ALGEZARES","ALJUCER",..: 16 11 7 7 33 14 12 34 34 34 ...
##  $ Barrio              : Factor w/ 72 levels "ALGEZARES","ALJUCER",..: 18 12 8 8 42 16 13 21 21 21 ...
##  $ Contador            : chr  "D12FE134047Z" "D12UF117329F" "D12UF117321X" "12-544140" ...
##  $ FechaInstContador   : Date, format: "2013-05-13" "2015-05-20" ...
##  $ Telectura           : Factor w/ 2 levels "N","S": 2 2 2 1 2 2 2 1 2 1 ...
##  $ TipoContador        : Factor w/ 5 levels "ELECTRONICO",..: 2 2 2 2 2 2 2 2 4 2 ...
##  $ ModeloContador      : Factor w/ 11 levels "ALTAIR V4 C",..: 7 7 7 4 7 7 7 7 5 4 ...
##  $ Calibre             : Factor w/ 11 levels "0","13","15",..: 7 8 8 3 8 7 7 7 5 6 ...
##  - attr(*, ".internal.selfref")=<externalptr>

3.2.2.2 Fichero 1.1: primera versión CONTRATOS_AGUA_PARQUES_Y_JARDINES.csv

paths$geo("contratos")
## [1] "./datosTFG/geo/CONTRATOS_AGUA_PARQUES_Y_JARDINES.csv"
contratos <- DatGeo$leerContratos(contratos = relaciones$Contrato)
  • Este fichero contiene 330 contratos de abastecimiento de agua con EMUASA por parte del ayuntamiento de Murcia, principalmente.
  • Cada contrato tiene un identificador único (Contrato), que es clave primaria.
  • El contrato depende del tipo de instalación (parque, jardín, …) sobre la que se hace (TipoContrato).

Vemos que están todos los contratos de relaciones que nos interesan.

str(contratos)
## Classes 'data.table' and 'data.frame':   330 obs. of  2 variables:
##  $ Contrato    : chr  "6351306" "6351960" "6353617" "6353635" ...
##  $ TipoContrato: Factor w/ 12 levels "ALTA AYUNTAMIENTO FUENTE",..: 3 3 3 3 3 3 3 3 3 3 ...
##  - attr(*, ".internal.selfref")=<externalptr>

3.2.3 Tarea 1. Reuniendo todos los datos geográficos.

3.2.3.1 Tarea 1.1. Qué zona riega cada contador.

Juntamos todas las tablas anteriores en una. 

dg <- DatGeo$reunirTablas()
## [2022-03-11 12:50:36] Reuniendo zonas verdes, acometidas y contratos...
str(dg)
## Classes 'sf' and 'data.frame':   362 obs. of  31 variables:
##  $ Acometida           : chr  "2182698" "2182749" "2182731" "2182694" ...
##  $ ZonasVerdes         : int  1 4 7 2 1 6 1 6 5 1 ...
##  $ Poligonos           : int  1 4 7 2 1 5 1 6 5 1 ...
##  $ MasGrande           : num  9.34 200.18 514.74 151.94 99.77 ...
##  $ MasPeque            : num  9.34 143.65 16.95 44.78 99.77 ...
##  $ Area                : num  9.34 714.54 1697 196.72 99.77 ...
##  $ MasLargo            : num  15.7 76.1 88.6 165.2 63.9 ...
##  $ MasCorto            : num  15.7 55 19.3 31.5 63.9 ...
##  $ Perimetro           : num  15.7 244.8 418.3 196.7 63.9 ...
##  $ MaxNodos            : int  16 18 19 26 8 30 10 9 46 7 ...
##  $ MinNodos            : int  16 14 5 9 8 17 10 5 5 7 ...
##  $ Nodos               : int  16 68 68 35 8 119 10 36 68 7 ...
##  $ Abcisa              : num  665925 665807 665384 665668 665551 ...
##  $ Norte               : num  4201430 4200964 4200542 4201059 4200921 ...
##  $ Contrato            : chr  "6393621" "6393645" "6393703" "6430007" ...
##  $ Actividad           : Factor w/ 4 levels "ACTIVIDAD GENERAL",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Ext/Int             : Factor w/ 2 levels "E","I": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Prop/Inqui          : Factor w/ 3 levels "C","I","P": 1 1 3 1 1 1 1 1 3 3 ...
##  $ Direccion           : chr  "PL.CANALEJAS" "CN.TERRERAS" "SANTUARIO DE LA FUENSANTA" "CALVARIO" ...
##  $ Indicacion          : Factor w/ 88 levels "10-JA","15 M-1-JA",..: 14 36 36 14 14 62 3 5 73 14 ...
##  $ ComplementoDireccion: chr  NA NA NA NA ...
##  $ Localidad           : Factor w/ 49 levels "ALGEZARES","ALJUCER",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Barrio              : Factor w/ 72 levels "ALGEZARES","ALJUCER",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Contador            : chr  "I18BB036286S" "D12UF128058K" NA "D11TA835359D" ...
##  $ FechaInstContador   : Date, format: "2018-10-05" "2013-01-15" ...
##  $ Telectura           : Factor w/ 2 levels "N","S": 2 2 1 2 1 2 2 2 2 1 ...
##  $ TipoContador        : Factor w/ 5 levels "ELECTRONICO",..: 4 2 NA 4 2 2 4 4 4 2 ...
##  $ ModeloContador      : Factor w/ 11 levels "ALTAIR V4 C",..: 5 7 NA 5 7 7 5 5 5 4 ...
##  $ Calibre             : Factor w/ 11 levels "0","13","15",..: 4 8 1 3 7 7 3 4 3 3 ...
##  $ TipoContrato        : Factor w/ 12 levels "ALTA AYUNTAMIENTO FUENTE",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ geometry            :sfc_MULTIPOLYGON of length 362; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:16, 1:2] 665925 665926 665927 665928 665928 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:30] "Acometida" "ZonasVerdes" "Poligonos" "MasGrande" ...

3.2.3.2 Tarea 1.2. Qué estación meteorológica está más cerca de cada jardín.

Añadimos la Distancia en \(m\) a la estación meteorológica más cercana (Principal) de la cual extraer las variables ambientales del parque. Además guardamos los identificadores de las otras estaciones para cada jardín, en orden por cercanía (AlternativaX, X=1,…,4).

# dg <- DatGeo$asignarEstMet(dg)
dg <- DatGeo$asignarEstMet(dg %>% dplyr::select())
## [2022-03-11 12:50:37] Asignando estaciones principal y alternativas a los parques...
str(dg, max.level = 1)
## Classes 'sf' and 'data.frame':   362 obs. of  7 variables:
##  $ geometry    :sfc_MULTIPOLYGON of length 362; first list element: List of 1
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  $ Distancia   : num  2114 1873 1511 1772 1658 ...
##  $ Principal   : chr  "MU62" "MU62" "MU62" "MU62" ...
##  $ Alternativa1: chr  "MU21" "MU21" "MU21" "MU21" ...
##  $ Alternativa2: chr  "MO12" "MO12" "MO12" "MO12" ...
##  $ Alternativa3: chr  "CA42" "CA42" "CA42" "CA42" ...
##  $ Alternativa4: chr  "CA91" "CA91" "CA91" "CA91" ...
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:6] "Distancia" "Principal" "Alternativa1" "Alternativa2" ...

3.2.3.3 Tarea 1.3. En qué mosaico está cada jardín.

Por último, obtenemos los identificadores de los mosaicos para Sentinel-2 en los que está contenido cada jardín (Mosaico, Moisaico2). Si sólo está en 1 ambos campos coinciden. Además contiene el mínimo número de mosaicos que hace falta para cubrirlo (MinMosaicos) y el número de mosaicos con los que interseca (MaxMosaicos). Vemos que:

  • todas las zonas están completamente contenidas (al menos) en un sólo mosaico;
  • la mayoría (337 / 363 = 92.84 %) están completamente contenidos en un sólo mosaico, mientras que unos pocos (26 / 363 = 7.16 %) están en la intersección de ambos mosaicos.
# system.time( dg <- DatGeo$asignarMosSen2(dg) )
system.time( dg <- DatGeo$asignarMosSen2(dg %>% dplyr::select()) )
## [2022-03-11 12:50:37] Identificando mosaicos Sentinel-2 en los que se encuentran los parques...
## Registered S3 method overwritten by 'geojsonsf':
##   method        from   
##   print.geojson geojson
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
##    user  system elapsed 
##   32.48    0.45   33.61
str(dg, max.level = 1)
## Classes 'sf' and 'data.frame':   362 obs. of  5 variables:
##  $ geometry   :sfc_MULTIPOLYGON of length 362; first list element: List of 1
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  $ MinMosaicos: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ MaxMosaicos: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Mosaico    : chr  "30SXH" "30SXH" "30SXH" "30SXH" ...
##  $ Mosaico2   : chr  "30SXH" "30SXH" "30SXH" "30SXH" ...
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA
##   ..- attr(*, "names")= chr [1:4] "MinMosaicos" "MaxMosaicos" "Mosaico" "Mosaico2"
dg %>% sf::st_drop_geometry() %>%
  dplyr::count(MinMosaicos, MaxMosaicos, Mosaico, Mosaico2, name = "Numero de parques")

3.2.3.4 Todos los pasos juntos (los métodos finales)

El método principal es DatGeo$initJardines(). Además de generar el fichero con los datos geográficos de todos los jardines juntos, los guarda en una carpeta, un fichero para cada parque, para la descarga de imágenes de Sentinel-2.

DatGeo$initJardines()
dg <- DatGeo$leerJardines()
str(dg, max.level = 1)
## Classes 'sf' and 'data.frame':   362 obs. of  41 variables:
##  $ Acometida           : chr  "2182698" "2182749" "2182731" "2182694" ...
##  $ ZonasVerdes         : int  1 4 7 2 1 6 1 6 5 1 ...
##  $ Poligonos           : int  1 4 7 2 1 5 1 6 5 1 ...
##  $ MasGrande           : num  9.34 200.18 514.74 151.94 99.77 ...
##  $ MasPeque            : num  9.34 143.65 16.95 44.78 99.77 ...
##  $ Area                : num  9.34 714.54 1697 196.72 99.77 ...
##  $ MasLargo            : num  15.7 76.1 88.6 165.2 63.9 ...
##  $ MasCorto            : num  15.7 55 19.3 31.5 63.9 ...
##  $ Perimetro           : num  15.7 244.8 418.3 196.7 63.9 ...
##  $ MaxNodos            : int  16 18 19 26 8 30 10 9 46 7 ...
##  $ MinNodos            : int  16 14 5 9 8 17 10 5 5 7 ...
##  $ Nodos               : int  16 68 68 35 8 119 10 36 68 7 ...
##  $ Abcisa              : num  665925 665807 665384 665668 665551 ...
##  $ Norte               : num  4201430 4200964 4200542 4201059 4200921 ...
##  $ Contrato            : chr  "6393621" "6393645" "6393703" "6430007" ...
##  $ Actividad           : Factor w/ 4 levels "ACTIVIDAD GENERAL",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Ext.Int             : Factor w/ 2 levels "E","I": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Prop.Inqui          : Factor w/ 3 levels "C","I","P": 1 1 3 1 1 1 1 1 3 3 ...
##  $ Direccion           : chr  "PL.CANALEJAS" "CN.TERRERAS" "SANTUARIO DE LA FUENSANTA" "CALVARIO" ...
##  $ Indicacion          : Factor w/ 88 levels "10-JA","15 M-1-JA",..: 14 36 36 14 14 62 3 5 73 14 ...
##  $ ComplementoDireccion: Factor w/ 36 levels "AGUA POTABLE PARA RIEGO",..: NA NA NA NA NA NA NA NA NA NA ...
##  $ Localidad           : Factor w/ 49 levels "ALGEZARES","ALJUCER",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Barrio              : Factor w/ 72 levels "ALGEZARES","ALJUCER",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Contador            : chr  "I18BB036286S" "D12UF128058K" NA "D11TA835359D" ...
##  $ FechaInstContador   : Date, format: "2018-10-05" "2013-01-15" ...
##  $ Telectura           : Factor w/ 2 levels "N","S": 2 2 1 2 1 2 2 2 2 1 ...
##  $ TipoContador        : Factor w/ 5 levels "ELECTRONICO",..: 4 2 NA 4 2 2 4 4 4 2 ...
##  $ ModeloContador      : Factor w/ 11 levels "ALTAIR V4 C",..: 5 7 NA 5 7 7 5 5 5 4 ...
##  $ Calibre             : Factor w/ 11 levels "0","100","13",..: 5 9 1 4 8 8 4 5 4 4 ...
##  $ TipoContrato        : Factor w/ 4 levels "ALTA AYUNTAMIENTO LOCALES",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Distancia           : num  2119 1877 1514 1776 1662 ...
##  $ Principal           : Factor w/ 5 levels "CA42","CA91",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ Alternativa1        : Factor w/ 5 levels "CA42","CA91",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Alternativa2        : Factor w/ 5 levels "CA42","CA91",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Alternativa3        : Factor w/ 5 levels "CA42","CA91",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Alternativa4        : Factor w/ 5 levels "CA42","CA91",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ MinMosaicos         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ MaxMosaicos         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Mosaico             : Factor w/ 2 levels "30SXG","30SXH": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Mosaico2            : Factor w/ 2 levels "30SXG","30SXH": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Geometry            :sfc_MULTIPOLYGON of length 362; first list element: List of 1
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "Geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:40] "Acometida" "ZonasVerdes" "Poligonos" "MasGrande" ...

3.2.4 Tarea 6.1.1. Dibujando simple features (sf)

Algunos ejemplos ilustrativos:

DatGeo$plotMapa(obj_sf=dg, group_var=NULL,
                color = "red", fillColor = "red")
DatGeo$plotMapa(obj_sf=dg, group_var=NULL,
                plot_var="Principal")
DatGeo$plotMapa(obj_sf=dg %>% dplyr::mutate(Distancia = Distancia / 1000),
                group_var=NULL,
                plot_var="Distancia",
                n = 5L)
DatGeo$plotMapa(obj_sf=dg, group_var="Barrio",
                filter_var="Localidad",
                filter_values=c("LA ALBERCA","SANTO ANGEL"),
                plot_var="Poligonos")
DatGeo$plotMapa(obj_sf=dg, group_var="Contrato",
                filter_var="Localidad",
                filter_values="LOBOSILLO",
                plot_var="Contrato",
                color = "black")
DatGeo$plotMapa(obj_sf=dg, group_var="Contrato",
                filter_var="Contrato",
                filter_values="6373577")
dgp <- sf::st_centroid(dg)
## Warning in st_centroid.sf(dg): st_centroid assumes attributes are constant over
## geometries of x
DatGeo$plotMapa(obj_sf=dgp, group_var=NULL)
DatGeo$plotMapa(obj_sf=dgp, group_var=NULL,
                cluster = F)
DatGeo$plotMapa(obj_sf=dgp, group_var=NULL,
                plot_var="Principal")
DatGeo$plotMapa(obj_sf=dgp %>% dplyr::mutate(Distancia = Distancia / 1000),
                group_var=NULL,
                plot_var="Distancia",
                n = 5,
                r = 5,
                cluster = F)
de <- DatGeo$leerEstMet() %>%
  sf::st_as_sf(coords = c("LONGITUD", "LATITUD"), crs = DatGeo$epsg$etrs89geo)
DatGeo$plotMapa(obj_sf=de, group_var=NULL)

Para guardar indicar un nombre_fichero:

# webshot::install_phantomjs()
m <- DatGeo$plotMapa(dg,"Contrato","Contrato","6373577",
                     nombre_fichero = "MapaCWorld")
## [2022-03-11 12:51:16] Guardando ./saves/en/MapaCWorld.png
m <- DatGeo$plotMapa(dg,"Contrato","Contrato","6373577",
                     nombre_fichero = "MapaCWorldCtrl", remove_controls = F)
## [2022-03-11 12:51:27] Guardando ./saves/en/MapaCWorldCtrl.png
m <- DatGeo$plotMapa(dg,"Contrato", width = 1050, height = 750,
                     nombre_fichero = "MapaWorld")
## [2022-03-11 12:51:32] Guardando ./saves/en/MapaWorld.png
m <- DatGeo$plotMapa(dg,"Contrato", width = 1050, height = 750,
                     nombre_fichero = "MapaWorldCtrl", remove_controls = F)
## [2022-03-11 12:51:44] Guardando ./saves/en/MapaWorldCtrl.png

3.2.5 Tarea 6.2. Visualización de datos (tablas e histogramas)

Para guardar indicar un f_save:

ply <- F
dflt$lan
## [1] "en"
# dflt$lan <- "sp"
f_save <- "png"
# f_save <- NULL
val <- 50; x_lab <- switch(dflt$lan,
  "sp" = "Número de polígonos",
  "en" = "Number of polygons"
)
DatInt$plotHist(dg, Poligonos, x_lab, sup = val,
                f_save = f_save, ply = ply)
## [2022-03-11 12:51:52] Eliminadas 2 observaciones de 362 (0.55%)
## [2022-03-11 12:51:52] Guardando Poligonos_l_50.png

# DatInt$plotHist(dg, Poligonos, x_lab, inf = val,
#                 f_save = f_save, ply = ply)
# gridExtra::grid.arrange(p1, p2, ncol=2)
DatInt$plotHist(dg, Poligonos, x_lab,
                f_save = f_save, ply = ply)
## [2022-03-11 12:51:52] Guardando Poligonos.png

n = dg %>% dplyr::filter(Poligonos == 1) %>% nrow()
N = dg %>% nrow()
msg(n," / ",N," = ",redondea2(n/N*100)," % (Poligonos = 1)")
## [2022-03-11 12:51:53] 111 / 362 = 30.66 % (Poligonos = 1)
val <- 12000; x_lab <- switch(dflt$lan,
  "sp" = latex2exp::TeX("Área total del parque ($m^2$)"),
  "en" = latex2exp::TeX("Total garden area ($m^2$)")
)
DatInt$plotHist(dg, Area, x_lab, sup = val,
                f_save = f_save, ply = ply)
## [2022-03-11 12:51:53] Eliminadas 5 observaciones de 362 (1.38%)
## [2022-03-11 12:51:53] Guardando Area_l_12000.png

DatInt$plotHist(dg, Area, x_lab,
                f_save = f_save, ply = ply)
## [2022-03-11 12:51:53] Guardando Area.png

val <- 10; x_lab <- switch(dflt$lan,
  "sp" = latex2exp::TeX("Distancia del parque a su estación principal ($km$)"),
  "en" = latex2exp::TeX("Distance from the garden to your main station ($km$)")
)
DatInt$plotHist(dg %>% dplyr::mutate(Distancia = Distancia/1000), Distancia,
                x_lab, f_save = f_save, ply = ply)
## [2022-03-11 12:51:54] Guardando Distancia.png

DatGeo$consultarDatGeo("tablaEstMet", dg = dg)
##              CA42 CA91 MO12 MU21 MU62 Total
## Principal      25    0   13   38  286   362
## Alternativa1    6   10   81  199   66   362
## Alternativa2   17   19  230   86   10   362
## Alternativa3  314   12   10   26    0   362
## Alternativa4    0  321   28   13    0   362
## Total         362  362  362  362  362  1810
escribirTablaLatex(DatGeo$consultarDatGeo("tablaEstMet", dg = dg), "estaciones")
## [2022-03-11 12:51:54] Guardando tabla ./saves/en/estaciones.txt
DatGeo$consultarDatGeo("tablaEstMet2", dg = dg)
##       CA42 CA91 MO12 MU21 MU62 Total
## CA42     0   10    0    0   15    25
## CA91     0    0    0    0    0     0
## MO12     0    0    0    0   13    13
## MU21     0    0    0    0   38    38
## MU62     6    0   81  199    0   286
## Total    6   10   81  199   66   362
escribirTablaLatex(DatGeo$consultarDatGeo("tablaEstMet2", dg = dg), "estaciones2")
## [2022-03-11 12:51:54] Guardando tabla ./saves/en/estaciones2.txt
t(DatGeo$consultarDatGeo("tablaMosSen2", dg = dg))
##       
##        30SXG 30SXH Ambos
##   [1,]    25   311    26
escribirTablaLatex(t(DatGeo$consultarDatGeo("tablaMosSen2", dg = dg)), "mosaicos")
## [2022-03-11 12:51:54] Guardando tabla ./saves/en/mosaicos.txt

4 Datos ambientales (DatAmb.R)

4.1 Fichero 4: Datos meteorológicos

4.1.1 Indicaciones

te paso el acceso a algunos datos del IMIDA que se llevan recogiendo desde 2015. Son 4 estaciones meteorológicas, cuya ubicación exacta se puede comprobar aquí: https://datosabiertos.regiondemurcia.es/carm/catalogo/medio-ambiente/estaciones-meteorologicas-de-la-red-siam-carm

  • MO12: Las Torres de Cotillas
  • MU62: Murcia
  • CA91: Fuente ÁLamo
  • CA42: Fuente Álamo.

El acceso: ftp://idei.imida.es/

Usario: um
Contraseña: ImidA_Siam2

Si fuera necesario o interesante creo que se podrían obtener los datos de cualquier estación.

Las variables son bastante intuitivas creo, te paso las que pueden no entenderse:

  • radmed: radiación media
  • radmax: radiación máxima
  • vvmed: velocidad del viento media
  • vvmax: velocidad del viento máxima
  • dvmed: dirección del viento media
  • prec: precipitaciones
  • dewpt: dew point (punto de rocio)
  • dpv: déficit de presión del vapor.

Por ahora no tenemos más estaciones meteorológicas, por lo que igual hay que coger 1 (o 2) para todos los parques. Digo 2 porque a veces habrán datos malos en una u otra.

4.1.2 Credenciales y valores por defecto

DatAmb$url
## [1] "ftp://idei.imida.es/"
DatAmb$userpwd
## [1] "um:ImidA_Siam2"

¡IMPORTANTE! Si se consiguen datos de más estaciones comprobar que realmente están accesibles (función DatAmb$consultarEstMet).

estaciones <- DatAmb$consultarEstMet()
estaciones
## [1] "CA42" "CA91" "MO12" "MU21" "MU62"
dflt$eemm
## [1] "CA42" "CA91" "MO12" "MU21" "MU62"

4.1.3 Descarga y preprocesamiento

El método principal es DatAmb$initDatMet(), que se apoya en DatAmb$downloadDatMet() y DatAmb$completarDatMet().

DatAmb$initDatMet()

¡IMPORTANTE! En DatAmb$completarDatMet se van completando las variables hasta que no faltan valores para DatAmb$vars, que son las que se utilizan para calcular la \(ET_0\).

Cuando actualizamos los datos realmente no se “actualizan”, se “machacan” los anteriores (init vs update). Además, se utiliza un único modelo por variable y pareja de estaciones. Todo ello puede resultar en un problema de reproducibilidad de resultados, ya que entre una actualización y la siguiente los modelos serán parecidos pero no idénticos, por lo que los valores por los que se sustituyan los NAs tampoco serán iguales.

NO SERIA MUY DIFICIL CAMBIAR ESTO:

  1. Para “actualizar”:
  • En DatAmb$completarDatMet al principio se leen los ficheros de cada estación. Cambiarlo para que el primero sea la “última versión” (utilizar DatAm$leerDatMet) añadiéndole los nuevos datos.
  • ¡Cambiar cómo se obtiene el nombre del fichero resultante! Observar variable k del bucle de DatAmb$completarDatMet (puede que las nuevas observaciones no tengan datos NA…)
  • Escribir DatAmb$updateDatMet a partir de DatAmb$initDatMet quitando la primera línea que borra todos los ficheros y la comprobación if (is.null(fichero)).
  1. Para conseguir reproducibilidad de resultados:
  • En la función na.replace de DatAmb$completarDatMet iterar sobre los valores a predecir (cuyos índices son idx_pred) obteniendo un modelo lineal a partir de sólo los datos anteriores (posicionces menores de índices de idx_train).

Realmente lo importante es hacer lo segundo, no lo primero. Sin embargo, aunque ahora no se tarde demasiado en completar los datos, si sólo hacemos el segundo punto sí que puede que vaya a empezar a ser un poco lento…

4.1.4 Leyendo los ficheros

Para leer los datos de una estación:

dm <- DatAmb$leerDatMet("MU21")
str(dm)
## Classes 'data.table' and 'data.frame':   49889 obs. of  15 variables:
##  $ Fecha   : POSIXct, format: "2015-11-25 00:00:00" "2015-11-25 01:00:00" ...
##  $ Tmed    : num  7.2 6.3 7.3 9 10 11.4 9.9 9.1 9 11 ...
##  $ Tmax    : num  7.6 7.7 7.7 9.5 10.5 12.8 11 9.8 9.7 11.8 ...
##  $ Tmin    : num  6.6 5.4 7.1 8 9.5 10.1 9.1 8.7 8.4 10.2 ...
##  $ HRmed   : num  65.5 71.1 64.8 57.8 54.2 48.4 54 58.1 61 55.1 ...
##  $ HRmax   : num  69 76.6 65.7 61.9 55.5 52.8 56.9 59.9 63 57.9 ...
##  $ HRmin   : num  63.9 64.4 63.5 55.5 52.5 43.7 50.8 55 58.8 52.5 ...
##  $ RADmed  : num  0 0 0 0 0 0 0 0.5 16.5 51.5 ...
##  $ RADmax  : num  0 0 0 0 0 0 0 2.6 34.9 75.7 ...
##  $ VVmed   : num  0.8 0.8 1.8 2.2 2.3 2.2 2.1 2.1 2.1 2.1 ...
##  $ VVmax   : num  2.5 2.6 3.2 4.1 7.1 7.1 4.5 3.7 3.5 3.9 ...
##  $ DVmed   : num  239 213 244 246 251 ...
##  $ Prec    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ PtoRocio: num  1.2 1.4 1.1 1.1 1.1 0.9 1 1.3 1.8 2.3 ...
##  $ DPV     : num  0.4 0.3 0.4 0.5 0.6 0.7 0.6 0.5 0.4 0.6 ...
##  - attr(*, ".internal.selfref")=<externalptr>

El parámetro estacion de DatAmb$leerDatMet puede ser identificador o un vector de ids, donde el primer elemento es la estacion principal y el resto las alternativas.

Para leer los datos de varias estaciones:

dms <- DatAmb$readDatMet(list("MU21", c("MU21", "MU62")))
str(dms, max.level = 1)
## List of 2
##  $ MU21     :Classes 'data.table' and 'data.frame':  49889 obs. of  15 variables:
##   ..- attr(*, ".internal.selfref")=<externalptr> 
##  $ MU21_MU62:Classes 'data.table' and 'data.frame':  49889 obs. of  15 variables:
##   ..- attr(*, ".internal.selfref")=<externalptr>
sapply(dms, function(dm) sum(!stats::complete.cases(dm)))
##      MU21 MU21_MU62 
##       565         2
sapply(dms, function(dm) sapply(dm, numNA))
##          MU21 MU21_MU62
## Fecha       0         0
## Tmed      517         0
## Tmax      516         0
## Tmin      523         0
## HRmed     518         0
## HRmax     518         0
## HRmin     524         0
## RADmed    551         0
## RADmax    551         0
## VVmed     512         0
## VVmax     513         0
## DVmed     512         2
## Prec      512         0
## PtoRocio  518         0
## DPV       518         0

El parámetro estaciones de DatAmb$readDatMet debe ser una lista del tipo estacion.

4.1.5 Visualizando los modelos de regresión lineal

Para guardar indicar un f_save:

f_save <- "svg"
f_save <- "png"
# f_save <- NULL
DatAmb$plotModelDatMet(princ = "CA42", alter = "CA91", vr = "HRmed",
                       fin = "2021-05-15", f_save = f_save)
## [2022-03-11 12:51:59] Guardando CA42_CA91_HRmed_l_2021-05-15.png
## [[1]]
## 
## Call:
## stats::lm(formula = y ~ x, data = df_train)
## 
## Coefficients:
## (Intercept)            x  
##      7.0648       0.9867  
## 
## 
## [[2]]

DatAmb$plotModelDatMet(princ = c("CA42","CA91"), alter = "MU62", vr = "HRmed",
                       fin = "2021-05-15", f_save = f_save)
## [2022-03-11 12:52:06] Guardando CA42_CA91_MU62_HRmed_l_2021-05-15.png
## [[1]]
## 
## Call:
## stats::lm(formula = y ~ x, data = df_train)
## 
## Coefficients:
## (Intercept)            x  
##     14.1417       0.9642  
## 
## 
## [[2]]

4.1.6 Tarea 6.2. Visualización de datos (tablas e histogramas)

DatAmb$consultarDatMet("tablaNA")

4.2 Paquete ET.PenmanMonteith

4.2.1 Valores por defecto

¡IMPORTANTE! La variable dflt$crops contiene los posibles valores de dflt$crop. En teoría da igual cuál elijamos, pues la diferencia es una constante. Comprobado experimentalmente para MU62.

dflt$crops
## [1] "short" "tall"
dflt$crop
## [1] "tall"

et0s

4.2.2 Preprocesamiento

Para instalar un paquete desde github:

devtools::install_github("VicenteYago/ET.PenmanMonteith")

El método principal es DatAmb$initET0(), que se apoya en DatAmb$calcularET0() que a su vez llama a ET.PenmanMonteith::et0().

DatAmb$initET0()

IMPORTANTE De nuevo, cuando actualizamos los datos realmente no se “actualizan”, se “machacan” los anteriores (init vs update). De nuevo, esto no sería muy difícil de cambiar… siempre y cuando ET.PenmanMonteith::et0() devuelva lo mismo (ver testET0.zip)

4.2.3 Leyendo uno de los ficheros

et0 <- DatAmb$leerET0(c("MU21","MU62"), "tall")
str(et0)
## Classes 'data.table' and 'data.frame':   2078 obs. of  2 variables:
##  $ Fecha: Date, format: "2015-11-25" "2015-11-26" ...
##  $ ET0  : num  3.25 3.89 1.63 2 1.2 ...
##  - attr(*, "param")=List of 1
##   ..$ crop: chr "tall"

4.3 Juntando todos los datos ambientales

Finalmente, los datos meteorológicos se agregan mediante el método estadístico correspondiente y se juntan con los valores de evapotranspiración de referencia (nada de esto se almacena, sino que se obtiene directamente).

da <- DatAmb$leerDatAmb("MU62", "tall")
str(da)
## Classes 'data.table' and 'data.frame':   2078 obs. of  16 variables:
##  $ Fecha   : Date, format: "2015-11-25" "2015-11-26" ...
##  $ Tmed    : num  14.3 17.9 14.5 12.6 11.6 ...
##  $ Tmax    : num  20.6 22 19.7 21.1 19.8 19.7 19.8 23 19.6 19.7 ...
##  $ Tmin    : num  4.8 14.1 8.9 6 5.4 4.2 4.5 5.3 6.3 5.7 ...
##  $ HRmed   : num  43.6 25.5 53.3 62.1 69.7 ...
##  $ HRmax   : num  72.6 42.4 80.9 87.4 89 91.3 90.1 90.9 89.7 91.2 ...
##  $ HRmin   : num  28.4 15.2 33.1 32.4 37.9 41.6 36 19.8 35.8 38.2 ...
##  $ RADmed  : num  122 133 128 128 129 ...
##  $ RADmax  : num  690 544 559 538 542 ...
##  $ VVmed   : num  1.233 1.692 0.758 0.692 0.537 ...
##  $ VVmax   : num  9.5 10.1 4.8 4.2 4.5 2.5 2 2.5 4.2 2.4 ...
##  $ DVmed   : num  247 271 164 211 210 ...
##  $ Prec    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ PtoRocio: num  2.01 -2.19 5.07 5.51 6.18 ...
##  $ DPVmed  : num  1.046 1.55 0.812 0.679 0.492 ...
##  $ ET0     : num  2.74 4.31 1.71 1.67 1.31 ...
##  - attr(*, ".internal.selfref")=<externalptr> 
##  - attr(*, "sorted")= chr "Fecha"
##  - attr(*, "param")=List of 1
##   ..$ crop: chr "tall"
# summary(da)

5 Datos de riego (Riego.R)

¡IMPORTANTE! El procesamiento de los datos de riego en run.R depende del servidor (dibulibu hace update y juno hace download). En vez de descargar los datos de dibulibu desde juno lo suyo sería pasar los ficheros de lecturas comprimidas y hacer todo en juno. Para ello, pasar los datos de la carpeta agua del servidor dibulibu al servidor juno en RiegoInteligente/datos/riego/agua, mediante rsync. Para cambiar el procesamiento de los datos de riego a juno cambiar Riego$nodename de "dibulibu" a "juno" y las credenciales Riego$host y Riego$pwd).

Hay días que no hay ficheros ZIP, pero se podrían obtener (ver CORREO Re: Posible entrevista con experto para dudas del TFG de Pedro J. Fernández Ruiz):

Riego$consultarFilesAgua("fechasFalt")
## [2022-03-11 12:52:18] Listando ficheros de lecturas comprimidas...
## [2022-03-11 12:52:19] Guardando listado de ficheros de lecturas comprimidas...

5.1 Fichero 3: Lecturas

5.1.1 Indicaciones

Todos los zip que tenemos se encuentran accesibles en un repositorio al que puedes acceder por ssh, yo hago en la terminal lo siguiente:

ssh aurora@155.54.204.34

y la contraseña es:

aauroralegustaelagua

Los ficheros se encentran en la carpeta “agua” y no te recomiendo que te los descargues todos, sino que cuando consigas simplifciar o agregar uno, hagamos una llamada a cada uno de esos zip pero sin descargarlos.

Esto sirve para ver qué fechas tenemos monitorizadas (para ver qué imágenes necesitamos y demás)

5.1.2 Credenciales

Riego$nodename
## [1] "dibulibu"
Riego$host
## [1] "aurora@155.54.204.34"
Riego$pwd
## [1] "aauroralegustaelagua"

5.1.3 Descarga o preprocesamiento

En el servidor dibulibu se ejecuta Riego$updateLecturas(), que se apoya en Riego$extraerLecturas().

Para ejecutar localmente basta llamar a Riego$downloadLecturas(), que es lo que hace juno.

Riego$downloadLecturas()

¡IMPORTANTE! Cuando actualizamos los datos de lecturas sí que se “actualizan” (gracias al fichero paths$proc()). Para “resetear” borrar dicho fichero y la carpeta paths$lects().

5.1.4 Leyendo uno de los ficheros

  • Los ficheros “lecturas/<contrato>.csv” contienen las lecturas del/de los contador/es del parque asociado a <contrato>.
  • El identificador del contador es Contador.
  • La fecha y hora a la que se ha hecho una telelectura es Fecha.
  • La telelectura del contador es Lectura. La unidad es el litro.
  • Adicionalmente se añaden variables que indican
    • si de una telelectura a la siguiente se produce un cambio de contador (CambiaContador),
    • el tiempo que pasa (Tiempo) y caúl sería el consumo (Consumo) entre una lectura y la siguiente.

Las muestras no están tomadas uniformemente. Tenemos datos desde el 29 de septiembre de 2017 hasta la actualidad, en general.

lectura <- Riego$leerLecturas("6397830", debug = T)
## [2022-03-11 12:52:19] 6397830: 34 lecturas erroneas
str(lectura)
## Classes 'data.table' and 'data.frame':   16808 obs. of  7 variables:
##  $ Expl2         : chr  "EMUASA-MURCIA" "EMUASA-MURCIA" "EMUASA-MURCIA" "EMUASA-MURCIA" ...
##  $ Contador      : chr  "00TC503427" "00TC503427" "00TC503427" "00TC503427" ...
##  $ Fecha         : POSIXct, format: "2017-09-29 00:38:29" "2017-09-29 01:38:35" ...
##  $ Lectura       : int  601309 601309 601309 601309 601309 601309 601309 601309 601309 601309 ...
##  $ Tiempo        :Formal class 'Duration' [package "lubridate"] with 1 slot
##   .. ..@ .Data: num  3606 3600 3600 3600 3598 ...
##  $ CambiaContador: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ Consumo       : int  0 0 0 0 0 0 0 0 0 0 ...

¡IMPORTANTE! a veces hay varias lecturas para un mismo instante. Si son iguales no pasa nada, se descartan indicándolo con un mensaje (debug = T). Cuando no son iguales parece que el contador se queda colgado y lo que hace es alternar entre dos valores durante varias tomas así que también se descartan y se notifica (debug = T). Con Riego$cleanLecturas() podemos eliminar de los ficheros las lecturas duplicadas (rmDup), las lecturas erróneas (rmErr) y quitar lecturas redundantes que no afectan al cálculo de los consumos (rmRed).

5.1.5 Algunas consultas (sobre resúmenes estadísticos)

5.1.5.1 Problema: Contratos que no tienen lecturas, tienen varios contadores o con contadores distintos en dg y en lecturas

Guardamos, para cada contrato, si no tiene, tiene uno, ninguno o varios contadores en lect. Los que no tienen ningún contador en son los que no tienen lecturas. Los que tienen varios es porque se lo han cambiado (no tiene varios simultáneamente). También guardamos, para los que tienen al menos uno en si es igual, casi igual o distinto al que hay en dg (si tiene varios nos fijamos en el último).

contadoresR <- Riego$leerResumenContadores()
# str(contadoresR)
Riego$consultarResumenContadores("sinLecturas", contadoresR) %>% nrow()
## [1] 47
Riego$consultarResumenContadores("variosContadores", contadoresR) %>% nrow()
## [1] 75

Es importante observar que a la hora de calcular los consumos tendremos que tener en cuenta cuándo se cambian los contadores.

5.1.5.2 ¿Con qué frecuencia y desde cuándo tenemos datos?

Para los contratos que tenemos lecturas hemos generado una tabla que guarda para cada contrato desde y hasta cuando tienen datos y un resumen sobre la frecuencia cuando hay consumo de las lecturas, porque será cuando realmente sesguemos los datos al calcular los consumos diarios.

lecturasR <- Riego$leerResumenLecturas()
# str(lecturasR)
# La mayoría (todos menos 62) de los contratos con lecturas recogen estos datos con telelectura desde el 29/09/2017
Riego$consultarResumenLecturas("tardias", "2017-09-30", res = lecturasR) %>% nrow()
## [1] 62
# Hay 12 contratos de los que empezamos a tener lecturas un poco tarde (después de `dflt$tardias`)
dflt$tardias
## [1] "2019-07-01"
Riego$consultarResumenLecturas("tardias", res = lecturasR) %>% nrow()
## [1] 12

La mayoría sigue recogiendo los datos hasta a día de hoy.

# Hay 17 contratos de los que no tenemos lecturas a día de hoy (hasta `dflt$noActuales`)
dflt$noActuales
## [1] "2021-05-08"
Riego$consultarResumenLecturas("noActuales", res = lecturasR) %>% nrow()
## [1] 17

Muchos de los contratos están como mucho 9 días sin telelectura en los periodos donde se registra un cambio de consumo.

# Hay 59 contratos de los que nos faltan muchas lecturas en un largo periodo de tiempo (`dflt$faltantesL`)
dflt$faltantesL
## [1] "9 days"
Riego$consultarResumenLecturas("faltantes", res = lecturasR) %>% nrow()
## [1] 59

5.1.5.3 Problema: Consumos negativos (reseteo, retorno o avería)

¿Cómo proceden cuando un contador alcanza la lectura máxima?

Cuando llega al máximo el contador mecánico empieza de nuevo por 0, en telelectura el contador no tiene el mismo límite por lo que se añade un dígito más, un contador de vueltas, los contadores que han dado la vuelta en facturación se ven con este dígito más que no tiene el contador mecánico, lo que no influye en la facturación ya que esta es a diferencias.

¿Por qué aparecen consumos de agua negativos, es decir, dos lecturas sucesivas de un mismo contador en las que la primera lectura es mayor que la siguiente?

Los consumos negativos aparecen por retornos en la red, por una avería en la válvula antirretorno el agua puede circular en poca cantidad en sentido inverso, debido a la diferencia de presiones. Pocos casos se han debido a avería del módulo de comunicaciones y se ha resuelto con el cambio del módulo. Otros casos de negativos se nos dan en contadores que se encuentran en los límites de sectores y el agua circula de forma habitual en ambos sentidos.

IMPORTANTE Cuando calculemos el consumo diario puede que muchas cantidades pequeñas negativas desaparezcan, pero no las grandes…

# Hay 126 contratos con consumos negativos
Riego$consultarResumenLecturas("negativo", res = lecturasR) %>% nrow()
## [1] 126

5.1.6 Tarea 6.2. Visualización de datos (tablas e histogramas)

contadoresR <- Riego$leerResumenContadores()
Riego$consultarResumenContadores("tabla", contadoresR)
##         <NA> Iguales CasiIguales Total
## Ninguno   47       0           0    47
## Uno        0      52         156   208
## Varios     0      74           1    75

Para guardar indicar un f_save:

ply <- F
dflt$lan
## [1] "en"
# dflt$lan <- "sp"
f_save <- "svg"
f_save <- "png"
# f_save <- NULL
lecturasR <- Riego$leerResumenLecturas() %>%
  dplyr::mutate_at(dplyr::vars(contains("Tiempo_")), lubridate::period_to_seconds) %>%
  dplyr::mutate_at(dplyr::vars(contains("Tiempo_")), as.difftime, units = "secs")
units(lecturasR$Tiempo_Min) <- "mins"
units(lecturasR$Tiempo_Mean) <- "hours"
units(lecturasR$Tiempo_Max) <- "days"
lecturasR %<>%
  dplyr::mutate_at(dplyr::vars(contains("Tiempo_")), as.double) %>%
  dplyr::mutate_at(dplyr::vars(contains("Tiempo_")), round, 2)
val <- as.Date("2017-09-30"); x_lab <- switch(dflt$lan,
  "sp" = "primera telelectura",
  "en" = "first remote reading"
)
DatInt$plotHist(lecturasR, Fecha_Min, x_lab, .bw = "month", inf = val,
                f_save = f_save, ply = ply)
## [2022-03-11 12:52:20] Eliminadas 221 observaciones de 283 (78.09%)
## [2022-03-11 12:52:20] Guardando Fecha_Min_Mes_ge_2017-09-30.png

DatInt$plotHist(lecturasR, Fecha_Min, x_lab, .bw = "quarter",
                f_save = f_save, ply = ply)
## [2022-03-11 12:52:20] Guardando Fecha_Min_Tri.png

val <- as.Date("2021-05-01"); x_lab <- switch(dflt$lan,
  "sp" = "última telelectura",
  "en" = "last remote reading"
)
DatInt$plotHist(lecturasR, Fecha_Max, x_lab, .bw = "month", sup = val,
                f_save = f_save, ply = ply)
## [2022-03-11 12:52:21] Eliminadas 266 observaciones de 283 (93.99%)
## [2022-03-11 12:52:21] Guardando Fecha_Max_Mes_l_2021-05-01.png

DatInt$plotHist(lecturasR, Fecha_Max, x_lab, .bw = "quarter",
                f_save = f_save, ply = ply)
## [2022-03-11 12:52:21] Guardando Fecha_Max_Tri.png

val <- 450; x_lab <- switch(dflt$lan,
  "sp" = "Longitud de intervalo de tiempo más largo con consumo pero sin lecturas (días)",
  "en" = "Length of longest time interval with consumption but without readings (days)"
)
DatInt$plotHist(lecturasR, Tiempo_Max, x_lab, .bw = "largo", sup = val,
                f_save = f_save, ply = ply)
## [2022-03-11 12:52:22] Eliminadas 15 observaciones NA de 283 (5.3%)
## [2022-03-11 12:52:22] Guardando Tiempo_Max_l_450.png

5.2 Consumos

5.2.1 Descarga o preprocesamiento

En el servidor dibulibu se ejecuta Riego$initConsumos(), que se apoya en Riego$calcularConsumos() y éste a su vez en Riego$deLecturasAConsumos()

Para ejecutar localmente podemos llamar a Riego$downloadConsumos(), que es lo que hace juno (o si hemos descargado las lecturas Riego$initConsumos()).

Riego$downloadConsumos()

¡IMPORTANTE! De nuevo, los consumos realmente no se “actualizan”, se “machacan” los anteriores (init vs update).

EN ESTE CASO CAMBIAR ESTE COMPORTAMIENTO PUEDE SER UN POCO MÁS COMPLICADO ya que hay que tener una lectura anterior y otra posterior al día, por lo que no podemos guardar el consumo de un día hasta que no tengamos una lectura del día siguiente… y todo depende del supuesto de que los datos de los ficheros ZIP vengan ordenados. No obstante, hay un intento de ello (Riego$updateConsumos() => Riego$calcularConsumos2() => Riego$deLecturasAConsumosPeriodo()).

¡IMPORTANTE! Hay dos versiones para calcular los consumos.

  • Versión 1: los consumos diarios entre cambios de contadores se establecen NA.
  • Versión 2 (por defecto): en este caso se calculan como se explica en el TFG.

5.2.2 Leyendo uno de los ficheros

  • Los ficheros “consumos/<contrato>.csv” contienen los consumos DIARIOS totales (ConsumoTotal) y por \(m^2\) (Consumo) del parque asociado a <contrato>.
  • La unidad es la misma que las lecturas (\(L\) o \(L/m^2\)).
consumo <- Riego$leerConsumos("6351960")
str(consumo)
## Classes 'data.table' and 'data.frame':   1395 obs. of  3 variables:
##  $ Fecha       : Date, format: "2017-09-30" "2017-10-01" ...
##  $ ConsumoTotal: num  332 460 5180 90 5200 ...
##  $ Consumo     : num  0.153 0.2122 2.3898 0.0415 2.399 ...

5.2.3 Algunas consultas (sobre resúmenes estadísticos)

5.2.3.1 Problema: Cambio de contadores (version 1) y consumos negativos

consumosR <- Riego$leerResumenConsumos()
# str(consumosR)
# Hay 52 contratos (de los 75 que cambian de contador) que nos faltan datos de consumo
Riego$consultarResumenConsumos("faltantes", 0, res = consumosR) %>% nrow()
## [1] 52
# Hay 30 contratos con consumos negativos
Riego$consultarResumenConsumos("negativo", res = consumosR) %>% nrow()
## [1] 30

5.2.3.2 Observación: Contratos de los que tenemos lecturas, pero su consumo es 0

# Hay 10 contratos con consumo 0
Riego$consultarResumenConsumos("cero", res = consumosR) %>% nrow()
## [1] 10

6 Datos de verdor (Verdor.R)

6.1 Ímagenes Satelitales y de Índices de Verdor

6.1.1 Indicaciones

Cuando el area del poligono es tan fina que no cubre el espacio del pixel de 10x10mts hace cosas raras, pero es algo q no afecta mucho al resultado final.

Acceso:

ssh gardens@juno.inf.um.es

la contraseña:

riego_inteligente

6.1.2 Credenciales y valores por defecto

Verdor$nodename
## [1] "juno"
Verdor$host
## [1] "gardens@juno.inf.um.es"
Verdor$pwd
## [1] "riego_inteligente"

¡IMPORTANTE! La variable dflt$inds contiene ALGUNOS posibles valores de dflt$ind. (con sen2r::list_indices(all = T) podemos ver todos los índices disponibles). Análogamente, dflt$masks contiene ALGUNOS posibles valores de dflt$mask, que incluyen todas las máscaras predefinidas por sen2r y un ejemplo (el último) de cómo definir una máscara a partir de las posibles etiquetas (números) del mapa SCL. El valor NA es especial, e indica sin aplicar máscara.

Doc ESA generación del SCL y significado: https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/level-2a/algorithm

Doc sen2r máscaras disponibles por defecto: https://sen2r.ranghetti.info/reference/s2_mask.html

dflt$inds
## [1] "NDVI" "SAVI" "EVI"  "ARVI" "NDMI"
dflt$ind
## [1] "NDVI"
dflt$masks
## [1] NA                   "nomask"             "nodata"            
## [4] "cloud_high_proba"   "cloud_medium_proba" "cloud_and_shadow"  
## [7] "clear_sky"          "land"               "scl_7_8_9"
dflt$mask
## [1] "scl_7_8_9"

6.1.3 Descarga o preprocesamiento en el servidor (Sección importante)

6.1.3.1 Descargando imágenes de Sentinel-2

  • Es necesaria una cuenta (hay 3, por lo que se pueden usar hasta 3 instancias)
  • El Espacio de trabajo es paths$sen2r() y ahí dentro cada cuenta tiene su propio espacio que incluye un fichero con las credenciales y donde se guardan los ficheros con los logs.
  • Es necesario guardar cada jardín en un GeoJSON diferente (**ya lo hace DatGeo$initJardines()**).
  • Verdor$updateImgInd() se encarga de repartir las fechas (que hay entre from y to) y las procesa en paralelo (n_accounts).
  • Cada instancia descarga sus imágenes (Verdor$downloadSentinel2()) usando las credenciales de una de las cuentas.
  • En realidad, sólo descarga las imágenes BOA y SCL, las imágenes de índices y su enmascaramiento se hace cuando se va a resumir la imagen correspondiente (más adelante).

6.1.3.2 ¿De qué jardines tenemos datos?

Éste es el listado que procesamos en el servidor (se trata de los jardines para los que conocemos su contrato y tenemos lecturas ACTUALES):

jardines <- dflt$jardinesSen2() %>% dplyr::select(Contrato)
str(jardines, max.level = 1)
## Classes 'sf' and 'data.frame':   266 obs. of  2 variables:
##  $ Contrato: chr  "6351960" "6353617" "6354706" "6355390" ...
##  $ Geometry:sfc_MULTIPOLYGON of length 266; first list element: List of 12
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "Geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
##   ..- attr(*, "names")= chr "Contrato"
# Todos los 266 jardines tienen carpeta:
Verdor$consultarDirsImgs("sinCarpeta") %>% nrow()
## [1] 0
nrow(jardines) == nrow(Verdor$consultarDirsImgs("conCarpeta"))
## [1] TRUE
# Sin embargo, 11 jardines no tienen nada dentro de su carpeta, por lo que tenemos imágenes de 255 contratos (pero ni idea de por qué pasa esto...)
Verdor$consultarDirsImgs("sinImagenes") %>% nrow()
## [1] 11

6.1.3.3 ¿Con qué frecuencia y desde cuando tenemos datos?

Tenemos datos cada 5 días (casi siempre) desde el 5 de julio de 2019 (Fecha_Min) para todos los contratos, lo que supone un máximo de 137 imágenes por contrato. El número de imágenes que faltan (Faltan) varían entre 4 y 5, dependiendo del jardín, pero en ningún caso falta más de 1 imagen seguida (Faltan_Max).

Verdor$consultarFilesImgsInd("total", ind = "BOA")
## [1] 174
imagesR <- Verdor$leerResumenFilesImgsInd(ind = "BOA")
summary(imagesR)
##    Contrato           Fecha_Min            Fecha_Max              Faltan     
##  Length:255         Min.   :2019-01-01   Min.   :2021-05-15   Min.   :4.000  
##  Class :character   1st Qu.:2019-01-01   1st Qu.:2021-05-15   1st Qu.:5.000  
##  Mode  :character   Median :2019-01-01   Median :2021-05-15   Median :5.000  
##                     Mean   :2019-01-01   Mean   :2021-05-15   Mean   :4.863  
##                     3rd Qu.:2019-01-01   3rd Qu.:2021-05-15   3rd Qu.:5.000  
##                     Max.   :2019-01-01   Max.   :2021-05-15   Max.   :5.000  
##    Faltan_Max
##  Min.   :1   
##  1st Qu.:1   
##  Median :1   
##  Mean   :1   
##  3rd Qu.:1   
##  Max.   :1

6.1.3.4 ¿Qué imágenes faltan? ¿Están disponibles online o en el LTA? ¿Cómo podemos solicitarlas? (SUPER IMPORTANTE)

Verdor$consultarFilesImgsInd("fechasFalt", fechasNo = NULL, ind = "BOA")
Verdor$consultarFilesImgsInd("tilesFalt", fechasNo = NULL, ind = "BOA",
                             fechas = as.Date("2021-03-31"))

Las 4 ó 5 imágenes que faltan son siempre las mismas: faltan las imágenes de 2019-05-21“,”2019-09-23“,”2020-10-22“, y”2021-05-10" para todos y las de “2021-03-31” para los que están en el tile “30SXH”.

Las causas de por qué faltan éstas fechas son o bien fallos o bien tareas de mantenimiento de Sentinel-2:

¡IMPORTANTE! En Verdor$fechasNoSen2 [ACTUALIZAR A MANO] guardamos estas fechas para obviarlas.

faltan <- Verdor$consultarFilesImgsInd("tilesFalt", ind = "BOA")
faltan
tl <- Verdor$checkAvailabilityDatesTiles(faltan)
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## Request failed [404]. Retrying in 1.4 seconds...
## Request failed [404]. Retrying in 1 seconds...
## Request failed [404]. Retrying in 1.5 seconds...
## Request failed [404]. Retrying in 1 seconds...
Verdor$cleanResponseAvailability(tl)

Para solicitar la disponibilidad necesitamos las credenciales de una de las cuentas:

so <- Verdor$orderLTA(tl)
## [2022-03-11 12:52:38] Check if products are already available for
##        download...
## [2022-03-11 12:52:38] Ordering 1 Sentinel-2 images stored in the Long
##        Term Archive...
## [2022-03-11 12:52:40] 1 of 1 Sentinel-2 images were correctly ordered.
##        You can check at a later time if the ordered products are
##        available online using the command:
##   safe_is_online("D:/Documentos/.sen2r/lta_orders/lta_20220311_125240.json")
str(so)
##  Named chr "https://apihub.copernicus.eu/apihub/odata/v1/Products('445d6c7c-7f99-4793-ac63-d0515c49d515')/$value"
##  - attr(*, "names")= chr "S2B_MSIL2A_20210331T104619_N0300_R051_T30SXH_20210331T134922.SAFE"
##  - attr(*, "available")= Named chr(0) 
##   ..- attr(*, "names")= chr(0) 
##  - attr(*, "notordered")= Named chr(0) 
##   ..- attr(*, "names")= chr(0) 
##  - attr(*, "path")= chr "D:/Documentos/.sen2r/lta_orders/lta_20220311_125240.json"
sen2r::safe_is_online(attr(so,"path"))
## 0 out of 1 products are online.
## S2B_MSIL2A_20210331T104619_N0300_R051_T30SXH_20210331T134922.SAFE 
##                                                             FALSE

Si no se han solicitado todas las imágenes, por ejemplo, porque hay un límite de solucitudes (creo que unas 20) por cuenta y hora, podemos solicitar después las que faltan, o usar otra de las cuentas:

Verdor$orderLTA(path = attr(so,"path"), account = 2)
## [2022-03-11 12:52:41] Check if products are already available for
##        download...
## [2022-03-11 12:52:41] Ordering 1 Sentinel-2 images stored in the Long
##        Term Archive...
## [2022-03-11 12:52:41] 1 of 1 Sentinel-2 images were correctly ordered.
##        You can check at a later time if the ordered products are
##        available online using the command:
##   safe_is_online("D:/Documentos/.sen2r/lta_orders/lta_20220311_125241.json")
##                                      S2B_MSIL2A_20210331T104619_N0300_R051_T30SXH_20210331T134922.SAFE 
## "https://apihub.copernicus.eu/apihub/odata/v1/Products('445d6c7c-7f99-4793-ac63-d0515c49d515')/$value" 
## attr(,"available")
## named character(0)
## attr(,"notordered")
## named character(0)
## attr(,"path")
## [1] "D:/Documentos/.sen2r/lta_orders/lta_20220311_125241.json"

IMPORTANTE Una vez online, las imágenes están unos 3 días disponibles, y la descarga y procesamiento de las imágenes es lento, así que tampoco interesa bombardear a solicitudes el servidor. Cuando yo estuve descargando imágenes (antes de actualizar el hardware de juno), lo normal era procesar unas 12 imágenes al día, es decir, se procesaban las imágenes de 1 mes (12 imágenes ~ 6 días de imágenes, ya que son 2 tiles ~ las imágenes de 1 mes, ya que la resolución temporal es de 5 días). Así que cada día solía hacer unas 12 peticiones. (Sí, estuve varias semanas, así que si queréis imágenes del LTA procurad automatizar esto, no hagáis lo que yo hice…)

6.1.4 Descarga o preprocesamiento local

Para ejecutar localmente primero ejecutar:

contrato <- '6373577'
Verdor$downloadImgInd(contrato,c("BOA","SCL"))

6.1.5 Leyendo una imagen

Para cargar una imagen TIFF como un RasterStack y así poder operar, dibujarla, etc.

contrato <- '6373577'
ind <- 'NDVI'; mask <- NA; dia <- '20210515'
img <- Verdor$leerImgInd(contrato, ind, mask, dia)
# str(img)

Para operar con sus valores lo normal es transforma a data.frame:

str( raster::as.data.frame(img) )
## 'data.frame':    228 obs. of  1 variable:
##  $ S2A2A_20210515_051_6373577_NDVI_10: num  NA NA NA NA NA NA NA NA NA NA ...
str( raster::as.data.frame(img, xy = TRUE) )
## 'data.frame':    228 obs. of  3 variables:
##  $ x                                 : num  664175 664185 664195 664205 664215 ...
##  $ y                                 : num  4205395 4205395 4205395 4205395 4205395 ...
##  $ S2A2A_20210515_051_6373577_NDVI_10: num  NA NA NA NA NA NA NA NA NA NA ...
summary( raster::as.data.frame(img, xy = TRUE) )
##        x                y           S2A2A_20210515_051_6373577_NDVI_10
##  Min.   :664175   Min.   :4205215   Min.   :0.1324                    
##  1st Qu.:664203   1st Qu.:4205255   1st Qu.:0.5337                    
##  Median :664230   Median :4205305   Median :0.7134                    
##  Mean   :664230   Mean   :4205305   Mean   :0.6722                    
##  3rd Qu.:664258   3rd Qu.:4205355   3rd Qu.:0.8403                    
##  Max.   :664285   Max.   :4205395   Max.   :0.8766                    
##                                     NA's   :163

6.1.6 Tarea 6.1.2. Visualizando archivos TIF (raster)

Para guardar indicar un f_save:

f_save <- "png"
f_save <- NULL

Podemos dibujar la imagen RGB, los índices o un histograma:

Verdor$plotImgInd(contrato, 'RGB', NA, dia, 0, dg, f_save = f_save)

Verdor$plotImgInd(contrato, 'SCL', NA, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

# Verdor$plotImgInd(contrato, dflt$ind, NA, dia, 1, dg, f_save = f_save)
# Verdor$plotImgInd(contrato, dflt$ind, NA, dia, 2, dg, f_save = f_save)
Verdor$plotImgInd(contrato, dflt$ind, NA, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

Verdor$plotImgInd(contrato, dflt$ind, "nomask", dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

Verdor$plotImgInd(contrato, dflt$ind, "nodata", dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

Verdor$plotImgInd(contrato, dflt$ind, dflt$mask, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

# Verdor$plotHistImgInd(contrato, dflt$ind, NA, dia, f_save = f_save)

6.1.7 Tenemos un problema con la detección de las nubes y las máscaras que van a producir outliers…

dia <- '20200902' # Nube blanca clasificada como suelo (amarillo)
Verdor$plotImgInd(contrato, 'RGB', NA, dia, 0, dg, f_save = f_save)

Verdor$plotImgInd(contrato, 'SCL', NA, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

dia <- '20210110' # Nube blanca clasificada como nieve (rosa) y zona oscura
Verdor$plotImgInd(contrato, 'RGB', NA, dia, 0, dg, f_save = f_save)

Verdor$plotImgInd(contrato, 'SCL', NA, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

dia <- '20191217' # Nube gris clasificada como cirro (azul)
Verdor$plotImgInd(contrato, 'RGB', NA, dia, 0, dg, f_save = f_save)

Verdor$plotImgInd(contrato, 'SCL', NA, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

dia <- '20200410' # ¿Nube gris o vegetacion?
Verdor$plotImgInd(contrato, 'RGB', NA, dia, 0, dg, f_save = f_save)

Verdor$plotImgInd(contrato, 'SCL', NA, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

dia <- '20200510'
Verdor$plotImgInd(contrato, 'RGB', NA, dia, 0, dg, f_save = f_save)

Verdor$plotImgInd(contrato, 'SCL', NA, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

dia <- '20200609'
Verdor$plotImgInd(contrato, 'RGB', NA, dia, 0, dg, f_save = f_save)

Verdor$plotImgInd(contrato, 'SCL', NA, dia, 3, dg, f_save = f_save)
## Warning: Use of `obj_df[[3]]` is discouraged. Use `.data[[3]]` instead.

6.1.8 Liberando espacio

Es preferible solicitar el nivel L2A al LTA que descargar un nivel L1C y transformarlo con sen2cor, por lo que no lo usamos y podríamos desinstalarlo…

s <- Verdor$conectarServVerdor()
status <- ssh::ssh_exec_wait(s, "./RiegoInteligente/disco.sh")
## gardens@juno:~$ df -h
## S.ficheros                    Tamaño Usados  Disp Uso% Montado en
## udev                            3,9G      0  3,9G   0% /dev
## tmpfs                           799M    81M  718M  11% /run
## /dev/mapper/Ubuntu16--vg-root    78G    31G   44G  42% /
## tmpfs                           3,9G      0  3,9G   0% /dev/shm
## tmpfs                           5,0M      0  5,0M   0% /run/lock
## tmpfs                           3,9G      0  3,9G   0% /sys/fs/cgroup
## /dev/vda1                       472M   212M  236M  48% /boot
## tmpfs                           799M      0  799M   0% /run/user/1001
## gardens@juno:~$ du -hs /home/gardens
## 4,4G /home/gardens
## gardens@juno:~$ du -hs /home/gardens/R
## 1,1G /home/gardens/R
## gardens@juno:~$ du -hs /home/gardens/RiegoInteligente
## 2,7G /home/gardens/RiegoInteligente
## gardens@juno:~$ du -hs /home/gardens/RiegoInteligente/code
## 12M  /home/gardens/RiegoInteligente/code
## gardens@juno:~$ du -hs /home/gardens/RiegoInteligente/code/sen2r
## 11M  /home/gardens/RiegoInteligente/code/sen2r
## gardens@juno:~$ du -hs /home/gardens/RiegoInteligente/datos/verdor/images
## 2,2G /home/gardens/RiegoInteligente/datos/verdor/images
## gardens@juno:~$ du -hs /home/gardens/RiegoInteligente/datos/verdor/summaries
## 37M  /home/gardens/RiegoInteligente/datos/verdor/summaries
status <- ssh::ssh_exec_wait(s, "du -hs .sen2r/sen2cor/")
## 596M .sen2r/sen2cor/
ssh::ssh_disconnect(s)

IMPORTANTISIMO!!!! NUNCA BORRAR LAS IMAGENES BOA NI SCL!! Tampoco es recomendable borrar las imágenes de índices (NDVI) sin enmascarar (NA), ya que si decidimos usar otra máscara, tendremos que volver a calcularlas. Se puede usar Verdor$cleanImgInd().

6.2 Resúmenes de imágenes de índices

6.2.1 Valores por defecto

¡IMPORTANTE! Verdor$funs contiene los valores que se calculan en Verdor$resumirImgInd(), pero modificar esta variable NO cambia los cálculos. Más adelante, tendremos que elegir uno de estos valores para representar los índices como serie unidimensional: dflt$agg, cuyos posibles valores son dflt$aggs. Debe ser un valor representativo/central; sin embargo, podría ser cualquier valor de Verdor$funs.

Verdor$funs
## [1] "Min"    "Q1"     "Median" "Mean"   "Q3"     "Max"    "n"
dflt$aggs
## [1] "Mean"   "Median"
dflt$agg
## [1] "Mean"

6.2.2 Descarga o preprocesamiento

Para ejecutar localmente primero llamar a Verdor$downloadSumInd().

En el servidor juno se ejecuta Verdor$updateSumInd(), que se apoya en Verdor$calcularSumInd(), que a su vez utiliza Verdor$resumirImgInd() y Verdor$getPathsImgInd(), que es la encargada de calcular las imágenes de índices y enmascararlas.

# Verdor$downloadSumInd('NDVI',c(NA,'scl_7_8_9'),contrato)
Verdor$downloadSumInd('NDVI',c(NA,'scl_7_8_9'))

En este caso los datos sí se “actualizan” (init vs update). Para “resetear” bastaría borrar la carpeta paths$sums(). Pero tened en cuenta que generar y resumir las imágenes tarda MUCHO tiempo.

6.2.3 Leyendo uno de los ficheros

idv <- Verdor$leerSumInd('6441469', 'NDVI', 'scl_7_8_9')
str(idv)
## Classes 'data.table' and 'data.frame':   144 obs. of  8 variables:
##  $ Fecha                : Date, format: "2019-07-05" "2019-07-10" ...
##  $ NDVI.scl_7_8_9.Min   : num  0.209 0.206 0.223 0.167 0.221 ...
##  $ NDVI.scl_7_8_9.Q1    : num  0.212 0.232 0.255 0.185 0.239 ...
##  $ NDVI.scl_7_8_9.Median: num  0.283 0.321 0.298 0.248 0.315 ...
##  $ NDVI.scl_7_8_9.Mean  : num  0.272 0.315 0.32 0.228 0.307 ...
##  $ NDVI.scl_7_8_9.Q3    : num  0.323 0.352 0.395 0.264 0.34 ...
##  $ NDVI.scl_7_8_9.Max   : num  0.352 0.442 0.455 0.299 0.439 ...
##  $ NDVI.scl_7_8_9.n     : int  9 9 9 9 9 9 9 9 9 9 ...
##  - attr(*, "param")=List of 2
##   ..$ ind : chr "NDVI"
##   ..$ mask: chr "scl_7_8_9"
  • Cada imagen se resume con 6 valores descriptivos (Min, Q1, Median, Mean,Q3, Max).
  • Además, guardamos el número de píxeles con datos (n).

6.2.4 Algunas consultas (sobre resúmenes estadísticos)

De los 255 contratos que tienen imágenes, hay 20 de los que en realidad no tenemos datos, es decir, que tienen imágenes “vacías” (que sólo tienen valores NA), por lo que sólo tenemos datos de 235 (pero ni idea de por qué pasa esto…).

# Hay 20 jardines que no tienen serie de índices sin enmascarar
Verdor$consultarResumenSumInd("sinSerie", mask = NA) %>% nrow()
## [1] 20

Al enmascarar perdemos los valores de otros 98-20=78 contratos; por lo que finalmente sólo tenemos datos de 157. IMPORTANTE Esto sí sabemos por qué pasa muchas veces… por la diferente resolución de las imágenes de índices (\(10\times10 m^2\)) y de los mapas SCL (\(20\times20 m^2\)). Podríamos no recortar los contornos para la adquisición de los mapas SCL…

# Hay 98 jardines que no tienen serie de índices enmascarada
Verdor$consultarResumenSumInd("sinSerie", mask = dflt$mask) %>% nrow()
## [1] 98

Realmente no debería ser necesario recortar los polígonos en los mapas SCL, bastaría recortarlos en las imgs BOA. Para ello habría que llamar 2 veces en vez de 1 a sen2r::sen2r() en Verdor$downloadSentinel2(). La primera, añadir el parámetro list_prods = "BOA" y la segunda los parámetros list_prods = "SCL", extent_as_mask = F.

Puesto que el recorte no es adecuado dejando fuera muchos pixels (es probable que incluso parques muy pequeños no recorten nada y sea esto lo que genera imágenes “vacías”), podría ser interesante que sen2r no recorte nada (simplemente añadir extent_as_mask = F sin las modificaciones anteriores) y hacer nosotros los recortes al resumir las imágenes (habría que modificar Verdor$resumirImgInd()).

Probar para algún caso y si merece la pena mover (no borrar) las BOA y/o los SCL para volver a descargar todas las imágenes…

# Hay 7 jardines que no parecen tener vegetacion (su NDVI no supera el valor `dflt$vegetacion`)
dflt$vegetacion
## [1] 0.3
Verdor$consultarResumenSumInd("sinVegetacion", mask = dflt$mask) %>% nrow()
## [1] 7

6.2.5 Tarea 6.2. Visualización de datos (tablas e histogramas)

Para guardar indicar un f_save:

ply <- F
dflt$lan
## [1] "en"
# dflt$lan <- "sp"
f_save <- "svg"
f_save <- "png"
# f_save <- NULL
idvsR <- Verdor$leerResumenSumInd() %>%
  dplyr::rename(N_NDVI = NDVI.scl_7_8_9.Mean_TotalNoNA) %>%
  dplyr::filter(N_NDVI > 0)
x_lab <- switch(dflt$lan,
  "sp" = "Valores NDVI disponibles",
  "en" = "Available NDVI values"
)
DatInt$plotHist(idvsR, N_NDVI, x_lab, f_save = f_save, ply = ply)
## [2022-03-11 12:53:21] Guardando N_NDVI.png

7 Tarea 5. Integración de las series de datos (DatInt.R)

7.1 Valores por defecto

¡IMPORTANTE! Se plantean dos maneras de integrar las series diarias con la serie de índices de verdor, idv, cuya frecuencia es de cinco días (dflt$ints):

  1. Redimensionar las series diarias (redim), es decir, de cada serie diaria se obtienen 5 series. Éste es el método por defecto (dflt$int).
  2. Agregar las series diarias (agg), es decir, cada 5 valores seguidos se agregan en 1.

En ambos casos se seleccionan las variables de interés (las del balance hídrico), se ajustan las fechas para que el valor idv de un día se relacione con los valores diarios de los 5 días anteriores (NO el actual) y se duplica “lagueada” la serie idv para relacionar un valor con el anterior, es decir, el de 5 días antes.

dflt$ints
## [1] "redim" "agg"
dflt$int
## [1] "redim"

De nuevo, en este caso los datos no se “actualizan” (init vs update)

7.2 Preprocesamiento

Para ejecutar localmente (si hemos ido descargando todos los datos) primero ejecutar DatInt$initDatInt(), que se apoya en DatInt$integrarSeries().

DatInt$initDatInt()

7.3 Leyendo los ficheros

Para leer los datos de uno de los contratos:

dts <- DatInt$leerDatInt('6441469', 'NDVI', 'scl_7_8_9', 'Mean', 'tall', 'redim')
str(dts)
## Classes 'data.table' and 'data.frame':   128 obs. of  18 variables:
##  $ Fecha                : Date, format: "2019-07-10" "2019-07-15" ...
##  $ Prec.5               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ET0.5                : num  7.49 7.5 8.37 7.59 8 ...
##  $ Consumo.5            : num  0 0 0 0.731 2.338 ...
##  $ Prec.4               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ET0.4                : num  6.99 7.5 7.5 8.43 7.29 ...
##  $ Consumo.4            : num  0 0 0 0.717 0.429 ...
##  $ Prec.3               : num  0 0 0 0 0 0 0 0.4 0 13.2 ...
##  $ ET0.3                : num  7.62 7.54 7.19 8.03 8.58 ...
##  $ Consumo.3            : num  0 0 0.85 1.324 0.422 ...
##  $ Prec.2               : num  0.1 0 0 0 0 0 0 0 0 0 ...
##  $ ET0.2                : num  4.99 7.41 7.52 8.25 7.28 ...
##  $ Consumo.2            : num  0 0 0.77 0.436 0.417 ...
##  $ Prec.1               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ET0.1                : num  8.36 8.09 8.17 8.29 6.94 ...
##  $ Consumo.1            : num  0 0 0.744 0.436 0.422 ...
##  $ NDVI.scl_7_8_9.Mean.5: num  0.272 0.315 0.32 0.228 0.307 ...
##  $ NDVI.scl_7_8_9.Mean  : num  0.315 0.32 0.228 0.307 0.292 ...
##  - attr(*, "param")=List of 5
##   ..$ ind : chr "NDVI"
##   ..$ mask: chr "scl_7_8_9"
##   ..$ agg : chr "Mean"
##   ..$ crop: chr "tall"
##   ..$ int : chr "redim"
dts <- DatInt$leerDatInt('6441469', 'NDVI', 'scl_7_8_9', 'Mean', 'tall', 'agg')
str(dts)
## Classes 'data.table' and 'data.frame':   128 obs. of  6 variables:
##  $ Fecha                  : Date, format: "2019-07-10" "2019-07-15" ...
##  $ Prec                   : num  0.1 0 0 0 0 0 0 0.4 0 13.2 ...
##  $ ET0                    : num  35.5 38 38.8 40.6 38.1 ...
##  $ Consumo                : num  0 0 2.36 3.64 4.03 ...
##  $ NDVI.scl_7_8_9.Mean    : num  0.272 0.315 0.32 0.228 0.307 ...
##  $ NDVI.scl_7_8_9.Mean.obj: num  0.315 0.32 0.228 0.307 0.292 ...
##  - attr(*, "param")=List of 5
##   ..$ ind : chr "NDVI"
##   ..$ mask: chr "scl_7_8_9"
##   ..$ agg : chr "Mean"
##   ..$ crop: chr "tall"
##   ..$ int : chr "agg"

Para leer los datos de todos los contratos (que tienen datos disponibles):

dts <- DatInt$readDatInt(NULL, 'NDVI', NA, 'Mean', 'tall', 'redim')
length(dts)
## [1] 235
dts <- DatInt$readDatInt(NULL, 'NDVI', 'scl_7_8_9', 'Mean', 'tall', 'redim')
length(dts)
## [1] 157

Tenemos una lista con nombres (los contratos) de data.frames. Para juntarlos en uno sólo añadiendo una columna Contrato:

str(dts[1:3], max.level = 1)
## List of 3
##  $ 6351960:Classes 'data.table' and 'data.frame':    128 obs. of  18 variables:
##   ..- attr(*, "param")=List of 5
##  $ 6353617:Classes 'data.table' and 'data.frame':    128 obs. of  18 variables:
##   ..- attr(*, "param")=List of 5
##  $ 6354706:Classes 'data.table' and 'data.frame':    128 obs. of  18 variables:
##   ..- attr(*, "param")=List of 5
str( dts <- listDF2DF(dts, "Contrato") )
## Classes 'data.table' and 'data.frame':   19768 obs. of  19 variables:
##  $ Contrato             : chr  "6351960" "6351960" "6351960" "6351960" ...
##  $ Fecha                : Date, format: "2019-07-10" "2019-07-15" ...
##  $ Prec.5               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ET0.5                : num  7.49 7.5 8.37 7.59 8 ...
##  $ Consumo.5            : num  4.554 4.507 4.673 0.157 4.77 ...
##  $ Prec.4               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ET0.4                : num  6.99 7.5 7.5 8.43 7.29 ...
##  $ Consumo.4            : num  0.286 4.517 4.641 0.125 4.812 ...
##  $ Prec.3               : num  0 0 0 0 0 0 0 0.4 0 13.2 ...
##  $ ET0.3                : num  7.62 7.54 7.19 8.03 8.58 ...
##  $ Consumo.3            : num  0.152 4.521 4.646 4.77 0.175 ...
##  $ Prec.2               : num  0.1 0 0 0 0 0 0 0 0 0 ...
##  $ ET0.2                : num  4.99 7.41 7.52 8.25 7.28 ...
##  $ Consumo.2            : num  4.54 0.12 7.949 5.808 0.129 ...
##  $ Prec.1               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ET0.1                : num  8.36 8.09 8.17 8.29 6.94 ...
##  $ Consumo.1            : num  4.503 0.166 4.664 4.84 6.302 ...
##  $ NDVI.scl_7_8_9.Mean.5: num  0.545 0.69 0.664 0.484 0.65 ...
##  $ NDVI.scl_7_8_9.Mean  : num  0.69 0.664 0.484 0.65 0.674 ...
##  - attr(*, ".internal.selfref")=<externalptr> 
##  - attr(*, "param")=List of 5
##   ..$ ind : chr "NDVI"
##   ..$ mask: chr "scl_7_8_9"
##   ..$ agg : chr "Mean"
##   ..$ crop: chr "tall"
##   ..$ int : chr "redim"

7.4 Algunas consultas (sobre resúmenes estadísticos)

7.4.1 Observaciones Totales, Completas y No Completas por contrato y entre todos

datosRna <- DatInt$leerResumenDatInt(mask = NA)
summary(datosRna)
##    Contrato           Obs_Total        Obs_CCs       Obs_NCCs     
##  Length:235         Min.   : 22.0   Min.   : 22   Min.   : 0.000  
##  Class :character   1st Qu.:128.0   1st Qu.:128   1st Qu.: 0.000  
##  Mode  :character   Median :128.0   Median :128   Median : 0.000  
##                     Mean   :126.4   Mean   :126   Mean   : 0.366  
##                     3rd Qu.:128.0   3rd Qu.:128   3rd Qu.: 0.000  
##                     Max.   :130.0   Max.   :130   Max.   :44.000
datosRna %>% dplyr::summarise_if(is.numeric, sum)
datosRscl <- DatInt$leerResumenDatInt()
summary(datosRscl)
##    Contrato           Obs_Total        Obs_CCs         Obs_NCCs    
##  Length:157         Min.   : 22.0   Min.   : 6.00   Min.   :  0.0  
##  Class :character   1st Qu.:128.0   1st Qu.:51.00   1st Qu.: 71.0  
##  Mode  :character   Median :128.0   Median :53.00   Median : 75.0  
##                     Mean   :125.9   Mean   :53.01   Mean   : 72.9  
##                     3rd Qu.:128.0   3rd Qu.:57.00   3rd Qu.: 77.0  
##                     Max.   :130.0   Max.   :69.00   Max.   :102.0
datosRscl %>% dplyr::summarise_if(is.numeric, sum)

7.5 Tarea 6.2. Visualización de datos (tablas e histogramas)

ply <- F
dflt$lan
## [1] "en"
# dflt$lan <- "sp"
f_save <- "svg"
f_save <- "png"
# f_save <- NULL
datosR <- DatInt$leerResumenDatInt() %>%
  dplyr::rename(N_Obs = Obs_CCs)
x_lab <- switch(dflt$lan,
  "sp" = "Observaciones completas disponibles",
  "en" = "Available complete cases"
)
DatInt$plotHist(datosR, N_Obs, x_lab, f_save = f_save, ply = ply)
## [2022-03-11 12:53:32] Guardando N_Obs.png

7.6 Tarea 6.3. Visualización de datos (gráficos de series temporales)

ply <- T
dflt$lan
## [1] "en"
# dflt$lan <- "sp"
f_save <- "svg"
f_save <- "png"
# f_save <- NULL
dts <- DatInt$leerSeries("6373577") %>%
  DatInt$reestructurarSeries()
str(dts)
## 'data.frame':    33428 obs. of  5 variables:
##  $ Fecha   : Date, format: "2015-11-25" "2015-11-26" ...
##  $ Entidad : chr  "MU62" "MU62" "MU62" "MU62" ...
##  $ Elemento: Factor w/ 8 levels "T","HR","RAD",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Variable: Factor w/ 22 levels "Tmed","Tmax",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Valor   : num  14.3 17.9 14.5 12.6 11.6 ...
dts %>% dplyr::distinct(Entidad)
dts %>% dplyr::distinct(Elemento, Variable) %>% dplyr::arrange(Elemento)
lapply(c("BH","ET0","Prec","Consumo","NDVI.scl_7_8_9.Mean"),
       function(v) DatInt$plotTS(dts, eleORvar = v, fin = dflt$fin,
                                 psz = 0.5, lsz = 0.3,
                                 f_save = f_save, ply = ply))
## [2022-03-11 12:53:33] Guardando BH_l_2021-05-17.png
## [2022-03-11 12:53:34] Guardando ET0_l_2021-05-17.png
## [2022-03-11 12:53:34] Guardando Prec_l_2021-05-17.png
## [2022-03-11 12:53:35] Guardando Consumo_l_2021-05-17.png
## [2022-03-11 12:53:35] Eliminadas 45 observaciones NA de 144 (31.25%)
## [2022-03-11 12:53:35] Guardando NDVI.scl_7_8_9.Mean_l_2021-05-17.png
## [[1]]
## 
## [[2]]
## 
## [[3]]
## 
## [[4]]
## 
## [[5]]